#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>

void poisson2d(int n, double **u, double xa, double ya, double h, double precision,
		double (* boundary_x1)(double x), double (* boundary_x2)(double x), 
		double (* boundary_y1)(double y), double (* boundary_y2)(double y), 
		double (* rho)(double x, double y) ){

	int i,j;
	double temp, change = precision +1;
	
	for(i=0;i<n;i++) u[i][0] = boundary_x1( xa + i*h );
	for(i=0;i<n;i++) u[i][n-1] = boundary_x2( xa + i*h );
	for(j=0;j<n;j++) u[0][j] = boundary_y1( ya + i*h );
	for(j=0;j<n;j++) u[n-1][j] = boundary_y2( ya + i*h );
	
	for(i=1;i<n-1;i++)
		for(j=1;j<n-1;j++)
			u[i][j]=0;
	u[n/2][n/2] = 1;
	
	while( change > precision ){
		
		change = 0;
		for(i=1;i<n-1;i++){
			for(j=1;j<n-1;j++){
				temp = u[i][j];
				if(i!=n/2 || j!=n/2)
				u[i][j] = 0.25 * (u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1] - h * h * rho(xa+i*h, ya+i*h) );
				change += fabs( u[i][j] - temp );
			}
		}
	}
	ofstream fp("poisson2d_point_charge.txt");
	for(j=0;j<n;j++){
		for(i=0;i<n;i++){
			fp<<u[i][j]<<"\t";
		}
		fp<<endl;
	}
	fp.close();
}

double x1(double x){
	return 0;
}
double x2(double x){
	return 0;
}
double y1(double y){
	return 0;
}
double y2(double y){
	return 0;
}
double rho(double x, double y){
	return 0;
}

int main(){

	int i,n=100;
	double xa=0, ya=0, h=1.0/n;

	double **u = new double *[n];
	for(i=0;i<n;i++) u[i] = new double [n];
	
	poisson2d( n, u, xa, ya, h, 0.001, x1, x2, y1, y2, rho);
	return 0;
}
